from __future__ import with_statement

import sys
import bisect
import math
import time
import itertools
import os.path
import logging
import optparse
import glob

import numpy
import numpy.random
import scipy.optimize
import scipy.stats

import common
import discrete_wlc
import refolding_common

logger = logging.getLogger("refolding")
logger.addHandler(logging.StreamHandler())
logger.setLevel(logging.INFO)

# Polyprotein SI: the real contour length of a folded (GB1) is ~32nm.
# folding rate constant is much more sensitive to force than the unfolding rate constant.
# Although GB1 folds really fast with a folding rate constant of ~700/s for GB1 at zero 
# force, the folding reaction of GB1 does not generate noticeable forces, as no force peaks 
# can be observed in the relaxation curve of polyprotein (GB1)8, in contrast with the ability 
# of ankyrin to generate force during folding.

slow_fallback = True


class Spec(object):

    def __init__(self, kf0, alpha):
        self.kf0 = kf0
        self.alpha = alpha

    def generate_time_to_first_folding(self, unfolded, f, size):
        common.check_gt(unfolded, 0)
        common.check_ge(f, 0)
        # min(X1, ..., Xn) is exp-distributed with lambda = sum(lambda_i)
        unfolding_rate = unfolded * self.kf(f)
        logger.debug("unfolding_rate=%r", unfolding_rate)
        if unfolding_rate == 0:
            # Avoid divistion by zero in computing the scale.
            t = numpy.repeat(numpy.inf, size)
        else:
            rv = scipy.stats.expon(loc=0, scale=1.0 / unfolding_rate)
            t = rv.rvs(size)
        return t
    
    def kf(self, f):
        return self.kf0 * numpy.exp(-self.alpha * f)

    def generate_folding_time_deltas(self, forces, size):
        """
        Row 0 is for waiting for one of the "unfolded" ones folding.
        Row 1 is for waiting for one of the "unfolded - 1" ones to fold.
        ...
        Row "unfolded - 1" is for waiting for the last one to fold.
        """
        common.check_nondecreasing(forces)
        deltas = numpy.empty((len(forces), size), dtype=refolding_common.float_type)
        for i in range(len(forces)):
            deltas[i, :] = self.generate_time_to_first_folding(len(forces) - i, forces[i], size)
        return deltas

    def __repr__(self):
        return "Spec(%r, %r)" % (self.kf0, self.alpha)

class Hover(object):

    def __init__(self, k, cl, z, unfolded):
        common.check_gt(cl, 0)
        common.check_ge(unfolded, 0)
        if z > cl:
            logger.warning("Very high extension: z=%r [nm] cl=%r [nm].", z * 1e9, cl * 1e9)
        self.k = k
        self.cl = cl
        self.z = z
        self.unfolded = unfolded

    def unfold(self, loop_length, modules_to_fold):
        return Hover(
            self.k, self.cl - modules_to_fold * loop_length, self.z, self.unfolded - modules_to_fold)

    def __repr__(self):
        return "Hover(k=%r, cl=%r, z=%r, unfolded=%r)" % (self.k, self.cl, self.z, self.unfolded)


class TimedHover(object):

    def __init__(self, hover, t, folded):
        common.check_ge(hover.unfolded, folded)
        self.hover = hover
        self.t = t
        self.folded = folded

    def unfolded(self):
        return self.hover.unfolded

    def z(self):
        return self.hover.z

    def unfolded_cl(self):
        return self.hover.cl

    def simplified_unfolded(self, loop_length):
        """
        From "Polyprotein of GB1 is an ideal artificial elastomeric protein":
        The total number of peaks in the first pulling curve between the lengths x and L0 is
        counted as Ntotal.
        """
        return min(self.unfolded(), (self.unfolded_cl() - self.z()) / loop_length)

    def __repr__(self):
        return "TimedHover(hover=%r, t=%r, folded=%r)" % (self.hover, self.t, self.folded)


class Model(object):

    def __init__(self, spec, kb_t_over_p, loop_length):
        """
        loop_length is zero for the inaccurate model
        """
        self.spec = spec
        self.kb_t_over_p = kb_t_over_p
        self.loop_length = loop_length

    def generate_folded_numbers(self, hover, t, size):
        forces = numpy.array([
                discrete_wlc.wlc_force(
                    self.kb_t_over_p,
                    hover.k * (hover.cl - i * self.loop_length),
                    hover.z / (hover.cl - i * self.loop_length))
                for i in range(hover.unfolded)])
        time_deltas = self.spec.generate_folding_time_deltas(forces, size)
        assert len(time_deltas) == hover.unfolded
        # We arbitrarily chose to consider one protein per column. Otherwise we could just as
        # well randomly permute each row before.
        folded = numpy.empty(size, dtype=int)
        for j in range(size):
            i, unfold_t = 0, 0.0
            while i < hover.unfolded:
                unfold_t += time_deltas[i, j]
                if unfold_t > t:
                    break
                i += 1
            folded[j] = i
        return folded

    def alpha(self):
        return self.spec.alpha

    def kf0(self):
        return self.spec.kf0

    def __repr__(self):
        return "Model(%r, %r, %r)" % (self.spec, self.kb_t_over_p, self.loop_length)

def _exp_int(e, fct, a, b):
    #return scipy.integrate.quad(lambda x : fct(x) * numpy.exp(e * x), a, b)[0]
    return scipy.integrate.romberg(lambda x : fct(x) * numpy.exp(e * x), a, b)

def slow_multiple_exp_integral(exponent, a, t):
    if len(exponent) == 0:
        return 1.0  # Should it be t - a ?
    e = exponent[1:]
    return _exp_int(
        exponent[0],
        lambda x : slow_multiple_exp_integral(e, x, t),
        a, t)

def _simplified_log_likelihood(unfolded, folded, t, kf):
    # Exact is False so that it works with numpy arrays.
    #logger.info("unfolded=%s folded=%s", list(unfolded), list(folded))
    combinations = scipy.misc.comb(unfolded, folded, exact=False)
    log_combinations = numpy.log(combinations)
    t_kf = t * kf
    log_folding_prob = folded * numpy.log(1 - numpy.exp(-t_kf))
    log_surviving_prob = -t_kf * (unfolded - folded)
    return log_combinations + log_folding_prob + log_surviving_prob

def compute_forces(kb_t_over_p, hovers):
    return numpy.array(
        [discrete_wlc.wlc_force(kb_t_over_p, hover.k * hover.cl, hover.z / hover.cl)
         for hover in hovers], dtype=refolding_common.float_type)

def _timed_hovers_info(model, timed_hovers):
    common.check_equal(model.loop_length, 0)
    hovers = [timed_hover.hover for timed_hover in timed_hovers]
    forces = compute_forces(model.kb_t_over_p, hovers)
    kf = model.spec.kf(forces)
    unfolded = numpy.array([hover.unfolded for hover in hovers])
    folded = numpy.array([timed_hover.folded for timed_hover in timed_hovers])
    t = numpy.array([timed_hover.t for timed_hover in timed_hovers])
    return hovers, forces, kf, unfolded, folded, t

def simplified_log_likelihood_from_timed_hovers(model, timed_hovers):
    hovers, forces, kf, unfolded, folded, t = _timed_hovers_info(model, timed_hovers)
    return simplified_log_likelihood(unfolded, folded, t, kf)

def simplified_log_likelihood(unfolded, folded, t, kf):
    slls = _simplified_log_likelihood(unfolded, folded, t, kf)
    return slls.sum()

def _simplified_partial_log_likelihood_partial_kf(unfolded, folded, t, kf):
    return t * (folded / (1.0 - numpy.exp(-t * kf)) - unfolded)

def simplified_partial_log_likelihood_partial_kf0(unfolded, folded, t, forces, kf0, alpha):
    exp_minus_alpha_forces = numpy.exp(-alpha * forces)
    kf = kf0 * exp_minus_alpha_forces
    return (exp_minus_alpha_forces * _simplified_partial_log_likelihood_partial_kf(
            unfolded, folded, t, kf)).sum()

def simplified_partial_log_likelihood_partial_alpha(unfolded, folded, t, forces, kf0, alpha):
    exp_minus_alpha_forces = numpy.exp(-alpha * forces)
    kf = kf0 * exp_minus_alpha_forces
    return (-kf0 * forces * exp_minus_alpha_forces * _simplified_partial_log_likelihood_partial_kf(
            unfolded, folded, t, kf)).sum()

def simplified_partial_2_log_likelihood_partial_kf0_2(unfolded, folded, t, forces, kf0, alpha):
    kf = kf0 * numpy.exp(-alpha * forces)
    exp_minus_t_kf = numpy.exp(- t * kf)
    return (-t ** 2 * folded * numpy.exp(-2 * alpha * forces) * exp_minus_t_kf / (1 - exp_minus_t_kf) ** 2).sum()

def simplified_estimate_kf0_from_events(kb_t_over_p, loop_length, events, alpha):
    return simplified_estimate_kf0_from_timed_hovers(
        kb_t_over_p, loop_length, [event.to_timed_hover() for event in events], alpha)

def estimate_kf0_from_events(kb_t_over_p, loop_length, events, alpha, kf0_slice):
    return estimate_kf0_from_timed_hovers(
        kb_t_over_p, loop_length, [event.to_timed_hover() for event in events], alpha, kf0_slice)

def simplified_estimate_kf0_from_timed_hovers(kb_t_over_p, loop_length, timed_hovers, alpha):
    unfolded = numpy.array(
        [timed_hover.simplified_unfolded(loop_length) for timed_hover in timed_hovers], dtype=int)
    folded = numpy.array([event.folded for event in timed_hovers], dtype=int)
    common.check_True((folded <= unfolded).all())
    t = numpy.array(
        [timed_hover.t for timed_hover in timed_hovers], dtype=refolding_common.float_type)
    hovers = [event.hover for event in timed_hovers]
    forces = compute_forces(kb_t_over_p, hovers)
    return simplified_estimate_kf0(unfolded, folded, t, forces, alpha)

def simplified_estimate_kf0(unfolded, folded, t, forces, alpha):
    logger.debug("simplified_estimate_kf0(..., alpha=%f)", alpha)
    fct = lambda kf0 : -simplified_log_likelihood(
        unfolded, folded, t, kf0 * numpy.exp(-alpha * forces))
    fct_prime = lambda kf0 : -simplified_partial_log_likelihood_partial_kf0(
        unfolded, folded, t, forces, kf0, alpha)
    fct_second = lambda kf0: -simplified_partial_2_log_likelihood_partial_kf0_2(
        unfolded, folded, t, forces, kf0, alpha)
    # fmin seems to be the most reliable and fastest method.
    optimisation_result = scipy.optimize.fmin(
        func=fct, x0=100.0,
        #fprime=fct_prime,
        #fhess=fct_second,
        maxiter=10000,
        full_output=True,
        disp=False,
        )
    xopt, fopt = optimisation_result[0:2]
    warnflag = optimisation_result[-1]
    common.check_equal(warnflag, 0)
    assert len(xopt) == 1
    result = xopt[0], -fopt
    logger.debug("result=%s", result)
    return result

simplified_delta_x_range = (0.1e-9, 10e-9)
simplified_alpha_points = 100


def estimate_kf0_from_timed_hovers(kb_t_over_p, loop_length, timed_hovers, alpha, kf0_slice):
    fct = lambda kf0 : numpy.inf if kf0 < 0 else -log_likelihood_from_timed_hovers(
        Model(Spec(kf0, alpha), kb_t_over_p, loop_length), timed_hovers)
    # scipy.optimize.fmin_powell
    optimisation_result = scipy.optimize.brute(
        fct, [kf0_slice],
        Ns=None, full_output=True,
        finish=_my_finish
        #finish=None
    )
    optimal_kf0, fval = optimisation_result[0:2]
    # 0-d array
    optimal_kf0 = refolding_common.float_type(optimal_kf0)
    optimal_ll = -fval
    logger.debug("optimal_kf0=%f optimal_ll=%f", optimal_kf0, optimal_ll)
    return optimal_kf0, optimal_ll

def _uci(params, kb_t_over_p, loop_length, timed_hovers):
    level = logging.DEBUG
    if logger.isEnabledFor(level):
        logger.log(level, "full inner params=%r", list(params))
    return numpy.inf if min(params) < 0 else -log_likelihood_from_timed_hovers(
        Model(Spec(params[0], params[1]), kb_t_over_p, loop_length), timed_hovers)

def estimate_from_timed_hovers(kb_t_over_p, loop_length, timed_hovers, alpha_slice, kf0_slice):
    fct = lambda params: _uci(params, kb_t_over_p, loop_length, timed_hovers)
    #fct = lambda params : -log_likelihood_from_timed_hovers(
    #    Model(Spec(params[0], params[1]), kb_t_over_p, loop_length), timed_hovers)
    # scipy.optimize.fmin_powell
    optimisation_result = scipy.optimize.brute(
        fct, [kf0_slice, alpha_slice],
        Ns=None, full_output=True,
        finish=_my_finish
        #finish=None
        )
    #xopt, fval = optimisation_result[0:2]
    xopt, fval, grid, jout = optimisation_result
    # 0-d array
    optimal_kf0, optimal_alpha = xopt
    optimal_ll = -fval
    logger.debug(
        "optimal_kf0=%r optimal_alpha=%r (optimal_x=%r nm) optimal_ll=%r",
        optimal_kf0, optimal_alpha, optimal_alpha * refolding_common.kb_t * 1e9, optimal_ll)
    return optimal_alpha, optimal_kf0, optimal_ll, (grid, -jout)

def estimate_from_events(kb_t_over_p, loop_length, events, alpha_slice, kf0_slice):
    return estimate_from_timed_hovers(
        kb_t_over_p, loop_length, [event.to_timed_hover() for event in events],
        alpha_slice, kf0_slice)

def _my_finish(*args, **kwargs):
    logger.info("Launching final optimisation.")
    return scipy.optimize.fmin_powell(*args, **kwargs)

# def _brute_finish(*args, **kwargs):
#     kwargs = kwargs.copy()
#     kwargs["maxiter"] = kwargs["maxfun"] = 5000
#     return scipy.optimize.fmin_powell(*args, **kwargs)

def simplified_estimate_from_timed_hovers(kb_t_over_p, loop_length, timed_hovers):
    unfolded = numpy.array(
        [timed_hover.simplified_unfolded(loop_length) for timed_hover in timed_hovers], dtype=int)
    folded = numpy.array([event.folded for event in timed_hovers], dtype=int)
    # It can fail becase of the simplified_unfolded stuff.
    common.check_equal(
        [timed_hovers[i] for i in range(len(folded)) if folded[i] > unfolded[i]], [])
    t = numpy.array([event.t for event in timed_hovers], dtype=refolding_common.float_type)
    hovers = [event.hover for event in timed_hovers]
    forces = compute_forces(kb_t_over_p, hovers)
    return simplified_estimate(unfolded, folded, t, forces)

def simplified_estimate(unfolded, folded, t, forces):
    alpha_range = simplified_delta_x_range[0] / refolding_common.kb_t, simplified_delta_x_range[1] / refolding_common.kb_t
    fct = lambda alpha : numpy.inf if alpha < 0 else -simplified_estimate_kf0(unfolded, folded, t, forces, alpha)[1]
    # We could pass finish=None to evaluate strictly only the grid.
    optimisation_result = scipy.optimize.brute(
        fct, [alpha_range], Ns=simplified_alpha_points, full_output=True,
        finish=_my_finish)
    optimal_alpha, fval = optimisation_result[0:2]
    optimal_kf0, optimal_sll = simplified_estimate_kf0(unfolded, folded, t, forces, optimal_alpha)
    common.check_equal(-fval, optimal_sll)
    logger.debug(
        "optimal_alpha=%r (optimal_x=%r nm) optimal_kf0=%r optimal_sll=%r",
        optimal_alpha, optimal_alpha * refolding_common.kb_t * 1e9, optimal_kf0, optimal_sll)
    return optimal_alpha, optimal_kf0, optimal_sll

def simplified_estimate_from_events(kb_t_over_p, loop_length, events):
    return simplified_estimate_from_timed_hovers(
        kb_t_over_p, loop_length, [event.to_timed_hover() for event in events])


class Movement(object):

    def __init__(self, z, t, unfolded, unfolded_cl, caught):
        self.z = z
        self.t = t
        self.unfolded = unfolded
        self.unfolded_cl = unfolded_cl
        self.caught = caught

    def __repr__(self):
        return "Movement(%r, %r, %r, %r, %r, %r)" % (
            self.z, self.t, self.unfolded, self.unfolded_cl, self.caught)

class Event(object):

    def __init__(self, k, z, t, unfolded, unfolded_cl, folded, name):
        self.k = k
        self.z = z
        self.t = t
        self.unfolded = unfolded
        self.unfolded_cl = unfolded_cl
        self.folded = folded
        self.name = name

    def __repr__(self):
        return "Event(%r, %r, %r, %r, %r, %r, %r)" % (
            self.k, self.z, self.t, self.unfolded, self.unfolded_cl, self.folded, self.name)

    def to_hover(self):
        return Hover(self.k, self.unfolded_cl, self.z, self.unfolded)

    def to_timed_hover(self):
        return TimedHover(self.to_hover(), self.t, self.folded)

    def with_t(self, t):
        return Event(self.k, self.z, t, self.unfolded, self.unfolded_cl, self.folded, self.name)


def multiple_exp_integral(exponent, t):
    """
    n = len(t)
    Computes integral t_0=0..t, t_1=t_0..t, ..., t_{n-1}=t_{n-2}..t
    exp(sum(e_i * t_i)) dt_0..dt_{n-1}.
    """
    n = len(exponent)
    common.check_equal((n, ), exponent.shape)
    b = numpy.ones((n + 1, n + 1), dtype=refolding_common.float_type)
    at_least_one_negative = False
    for j in range(1, n + 1):
        if at_least_one_negative:
            break
        cs = exponent[j - 1:].cumsum()
        for i in range(j, n + 1):
            s = cs[i - j]
            # On some computers we don't always get equality.
            assert abs(s - exponent[j - 1:i].sum()) < 1e-10, (s, exponent[j - 1:i].sum())
            numerator = numpy.exp(s * t) * b[j - 1, j - 1] - b[i, j - 1]
            b[i, j] = numerator / s
            common.check_False(numpy.isnan(b[i, j]), exponent)
            if math.copysign(1, b[i, j]) == -1:
                logger.debug(
                    "i=%d j=%d s=%r b=\n%r exponent=%r numpy.exp(s * t)=%r b[j - 1, j - 1]=%r b[i, j - 1]=%r",
                    i, j, s, b, exponent, numpy.exp(s * t), b[j - 1, j - 1], b[i, j - 1])
                # Rounding errors can make b negative.
                #assert abs(numerator) < 10, (numerator, s, b[i, j])
                at_least_one_negative = True
                break
    if at_least_one_negative:
        if not slow_fallback:
            return 0
        slow = slow_multiple_exp_integral(exponent, 0, t)
        logger.info("Bad: slow=%r exponent=%r", slow, exponent)
        result = slow
    else:
        result = b[n, n]
    return result

def log_likelihood_from_timed_hover(model, timed_hover):
    common.check_not_equal(timed_hover.unfolded(), 0, timed_hover)
    hovers = [timed_hover.hover.unfold(model.loop_length, i) for i in range(timed_hover.folded + 1)]
    assert timed_hover.folded + 1 == len(hovers)
    forces = compute_forces(model.kb_t_over_p, hovers)
    common.check_False((forces < 0).any(), forces)
    kf = model.spec.kf(forces)
    common.check_False((kf < 0).any(), kf)
    if (kf[:-1] == 0).any():
        # Rounding errors made a transition impossible.
        logger.warning("model=%r Zero folding rate", model)
        p = -numpy.inf
    else:
        common.check_nonincreasing(kf)
        unfolded = numpy.array(
            [timed_hover.unfolded() - i for i in range(timed_hover.folded + 1)],
            dtype=refolding_common.float_type)
        assert len(unfolded) == len(hovers)
        assert unfolded[-1] == timed_hover.unfolded() - timed_hover.folded
        unfolded_times_kf = unfolded * kf
        exponent = (numpy.roll(unfolded_times_kf, -1) - unfolded_times_kf)[:-1]
        common.check_False((exponent > 0).any())
        p = numpy.log(unfolded[:-1]).sum()
        common.check_False(numpy.isnan(p))
        p += numpy.log(kf[:-1]).sum()
        common.check_False(numpy.isnan(p))
        p -= unfolded_times_kf[-1] * timed_hover.t
        common.check_False(numpy.isnan(p))
        mei = multiple_exp_integral(exponent, timed_hover.t)
        logger.debug("exponent=%r kf=%r mei=%r", exponent, kf, mei)
        p += numpy.log(mei)
        logger.debug("final p=%r", p)
        common.check_False(numpy.isnan(p))
    if slow_fallback:
        common.check_False(p == numpy.inf, (model, timed_hover))
    return p

def log_likelihood_from_timed_hovers(model, timed_hovers):
    s = sum(
        (log_likelihood_from_timed_hover(model, timed_hover) for timed_hover in timed_hovers),
        refolding_common.zero)
    if slow_fallback:
        common.check_False(numpy.isinf(s))
    return s

def more_less_estimation(
    loop_length,
    fixed_alpha,
    events,
    estimate_kf0_from_events_function, estimate_from_events_function,
    method, log_likelihood_image, text):

    if fixed_alpha:
        kf0_fixed_alpha, ll_fixed_alpha = estimate_kf0_from_events_function(
            refolding_common.kb_t_over_p, loop_length, events, fixed_alpha)
        logger.info("%s kf0_fixed_alpha=%f ll_fixed_alpha=%r", method, kf0_fixed_alpha, ll_fixed_alpha)

    alpha, kf0, ll, (grid, jout) = estimate_from_events_function(
        refolding_common.kb_t_over_p, loop_length, events)
    logger.info(
        "%s: alpha=%f (delta_x=%f nm) kf0=%f ll=%r",
        method, alpha, alpha * refolding_common.kb_t * 1e9, kf0, ll)

    if log_likelihood_image:
        #print(grid.shape, grid[0], grid[1])
        #print(jout.shape, jout[0])
        
        import matplotlib
        matplotlib.use("PS")
        import pylab
        import matplotlib.backends.backend_ps
        import mpl_toolkits.mplot3d
        pylab.rc("font", size=5)

        # definitions for the axes
        left, width = 0.1, 0.65
        bottom, height = 0.1, 0.65
        bottom_h = left_h = left + width + 0.02
        rect_scatter = [left, bottom, width, height]
        rect_histx = [left, bottom_h, width, 0.2]
        rect_histy = [left_h, bottom, 0.2, height]
        
        figure = matplotlib.pyplot.figure(1, figsize=(8, 8))
        ax = mpl_toolkits.mplot3d.Axes3D(figure)
        #ax.plot_wireframe(grid[0], grid[1] * refolding_common.kb_t * 1e9, jout)
        ravelled_x = numpy.ravel(grid[0])
        ravelled_y = numpy.ravel(grid[1] * refolding_common.kb_t * 1e9)
        ravelled_jout = numpy.ravel(jout)
        #print("ravelled_jout.shape=%s" % (ravelled_jout.shape, ))
        colours = refolding_common.gradient(ravelled_jout)
        #print("colours=%r" % colours)
        ax.scatter3D(ravelled_x, ravelled_y, ravelled_jout, color=colours)
        ax.scatter3D([kf0], [alpha * refolding_common.kb_t * 1e9], [ll], color="black")
        if text:
            for (xx, yy, zz, colour) in zip(ravelled_x, ravelled_y, ravelled_jout, colours):
                t = matplotlib.axes.Axes.text(ax, xx, yy, "%.2f" % zz, color=colour)
                mpl_toolkits.mplot3d.art3d.text_2d_to_3d(t, zz, zdir=(1, 0, 0))
        t = matplotlib.axes.Axes.text(
            ax, kf0, alpha * refolding_common.kb_t * 1e9, "%.2f" % ll, color="black")
        mpl_toolkits.mplot3d.art3d.text_2d_to_3d(t, ll, zdir=(1, 0, 0))
        
        ax.set_xlabel("kf0")
        ax.set_ylabel("delta_x")
        ax.set_zlabel("log_likelihood")

        ax.set_xlim3d(min(kf0, grid[0].min()), max(kf0, grid[0].max()))
        ax.set_ylim3d(
            min(alpha, grid[1].min()) * refolding_common.kb_t * 1e9,
            max(alpha, grid[1].max()) * refolding_common.kb_t * 1e9)
        ax.set_zlim3d(jout.min(), jout.max())

        canvas = matplotlib.backends.backend_ps.FigureCanvasPS(figure)
        canvas.print_figure(log_likelihood_image)


def read_events(root_dir, name):
    full_path = os.path.join(root_dir, name + ".full")
    entries = refolding_common.read_entries(full_path)
    events = []
    for entry in entries:
        common.check_True(isinstance(entry, tuple))
        k, hover_nm, t_ms, unfolded, unfolded_cl_nm, folded, filename = entry
        events.append(
            Event(k, hover_nm * 1e-9, t_ms * 1e-3, unfolded, unfolded_cl_nm * 1e-9, folded, filename))
    return events

def check_names_unique(events):
    all_names_list = [event.name for event in events]
    all_names_set = set(all_names_list)
    common.check_False([name for name in all_names_set if all_names_list.count(name) != 1])

def add_reading_options(parser):
    parser.add_option(
        "--root_dirs", help="dirs having under them .full files") 
    parser.add_option(
        "--ignore_starred", choices=("True", "true", "False", "false"),
        help="ignore marked curves", default="True")

def is_ignore_starred(options):
    return {"false" : False, "true" : True}.get(options.ignore_starred.lower())

def read_all_events(root_dirs, ignore_starred):
    all_events = []
    for root_dir in root_dirs.split(":"):
        if root_dir:
            for filepath in glob.iglob(os.path.join(root_dir, "*.full")):
                name = os.path.basename(os.path.splitext(filepath)[0])
                events = read_events(root_dir, name)
                events = [event for event in events if not event.name.endswith("*")] if ignore_starred else events
                print("%d good curves in %s" % (len(events), os.path.join(root_dir, name)))
                all_events.extend(events)
    return all_events

def main(params):
    parser = optparse.OptionParser()
    add_reading_options(parser)
    parser.add_option("--xf_nm", help="xf if already known")
    parser.add_option("--loop_nm", help="loop length")
    parser.add_option("--random_sampling_with_replacement", choices=("true", "True", "false", "False"))
    parser.add_option("--random_sampling_time_separators_ms")
    parser.add_option("--xf_start_nm", help="")
    parser.add_option("--xf_end_nm", help="")
    parser.add_option("--xf_step_nm", help="")
    parser.add_option("--kf0_start", help="")
    parser.add_option("--kf0_end", help="")
    parser.add_option("--kf0_step", help="")
    parser.add_option("--log_likelihood_image", help="optional eps filename")
    parser.add_option("--text", choices=("True", "true", "False", "false"))

    options, args = parser.parse_args(params)
    if args:
        parser.print_help()
        return -1
    print(options)


    ignore_starred  = is_ignore_starred(options)
    all_events = read_all_events(options.root_dirs, ignore_starred)
    check_names_unique(all_events)
    logger.info("Total: %d curves", len(all_events))

    random_sampling_with_replacement = {"false" : False, "true" : True}.get(
        options.random_sampling_with_replacement.lower())
    random_sampling_time_separators_ms = options.random_sampling_time_separators_ms
    seps_ms = map(float, random_sampling_time_separators_ms.split(":")) if random_sampling_time_separators_ms else []
    seps = sorted([x * 1e-3 for x in seps_ms])
    common.check_equal(len(seps), len(set(seps)))
    if random_sampling_with_replacement:
        seed = numpy.random.randint(0, 99999999)
        logger.info("seed=%d", seed)
        numpy.random.seed(seed)

        buckets = {}
        for i in range(len(all_events)):
            event = all_events[i]
            bucket = bisect.bisect_left(seps, event.t)
            bucket_indices = buckets.setdefault(bucket, [])
            bucket_indices.append(i)
        assert range(len(all_events)) == sorted(sum(buckets.values(), []))
        sampled_buckets = {}
        for (time_index, bucket) in buckets.items():
            sampled_buckets[time_index] = [
                bucket[r] for r in numpy.random.randint(0, len(bucket), len(bucket))]
        logger.info("sampled_buckets=%s", repr(sampled_buckets))
        indices = sum(sampled_buckets.values(), [])
        all_events = [all_events[index] for index in indices]

    loop_length = refolding_common.float_type(options.loop_nm) / 1e9
    xf_start = refolding_common.float_type(options.xf_start_nm) / 1e9
    xf_end = refolding_common.float_type(options.xf_end_nm) / 1e9
    xf_step = refolding_common.float_type(options.xf_step_nm) / 1e9
    kf0_start = refolding_common.float_type(options.kf0_start)
    kf0_end = refolding_common.float_type(options.kf0_end)
    kf0_step = refolding_common.float_type(options.kf0_step)

    if options.xf_nm:
        xf_m = refolding_common.float_type(options.xf_nm) / 1e9
        fixed_alpha = xf_m / refolding_common.kb_t
    else:
        fixed_alpha = None

    more_less_estimation(
        loop_length,
        fixed_alpha,
        all_events,
        simplified_estimate_kf0_from_events,
        lambda *args : simplified_estimate_from_events(*args) + ((None, None), ),
        "Simplified", None, False
        )

    alpha_start = xf_start / refolding_common.kb_t
    alpha_end = xf_end / refolding_common.kb_t
    alpha_step = xf_step / refolding_common.kb_t

    alpha_slice = slice(alpha_start, alpha_end, alpha_step)
    kf0_slice = slice(kf0_start, kf0_end, kf0_step)

    log_likelihood_image = options.log_likelihood_image
    text = {"false" : False, "true" : True}.get(options.text.lower())

    more_less_estimation(
        loop_length,
        fixed_alpha,
        all_events,
        lambda *args : estimate_kf0_from_events(*(args + (kf0_slice,))),
        lambda *args : estimate_from_events(*(args + (alpha_slice, kf0_slice))),
        "Full", log_likelihood_image, text)


if __name__ == "__main__":
    main(sys.argv[1:])
